﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;
using NCM_MSTest.Alg;

namespace NCM_MSTest
{
    /// <summary>
    /// 数值积分与微分
    /// I(f) =\int {f(x)dx} [a,b]
    /// </summary>
    internal class NumbericalIntegrationDifferentiation
    {
        /*
        等距节点求积公式
        矩形求积公式
            I   = \int{f(x)dx}[a,b]
                = \sum{Si} = \sum{hi*f(xi)}
            Si = f(xi)*(xi - x_{i-1})
         */

        /// <summary>
        /// 矩形求积公式(具有一阶精度)
        /// I=\int_a^b{f(x)}dx=\sum_{i=1}^{n}{hi*f(xi)}
        /// I=SUM( hi * f(xi) ) i=1,...,n
        /// 将区间[a,b]分成n个等间距的小区间[xi,xj]...
        /// 在任意区间[xi,xj]内任取一点xi,对应函数f(xi)
        /// 则矩形面积Si = hi * f(xi)=(xj - xi)*f(xi)
        /// 则积分 I = SUM(Si)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="step"></param>
        /// <param name="I"></param>
        public static void Integral_Rectangle(F f, double a, double b
            , double step
            , out double I)
        {
            I = 0;
            int n = (int)Math.Floor((b - a) / step);
            double xi = 0;
            for (int i = 0; i < n; i++)
            {
                xi = a + step * i;
                I += f(xi);
                //Console.WriteLine($"{i}\t {xi:f4}\t {f(xi):f5} \t {I:f5}");
            }
            I = I * step;
            I += f(b) * (b - xi);
            //Console.WriteLine($" \t {I}");
        }

        /// <summary>
        /// 梯形求积公式(具有二阶精度)
        /// I=h/2*(f(a)+f(b)+ 2*SUM(f(xi)) ) i=1,...,n-1
        /// 小直角梯形的高：hi = (b - a)/n
        /// 小直角梯形的面积：Si = (f(xi)+f(xj))/2 * hi
        /// I= SUM(Si)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="h"></param>
        /// <param name="I"></param>
        public static void Integral_Trapezoidal(F f, double a, double b
            , double h
            , out double I)
        {
            I = double.NaN;
            int n = (int)Math.Floor((b - a) / h) + 1;
            double s = f(a) / 2 + f(b) / 2;
            double x = a;
            Console.WriteLine($"i \t x \t f(x)");
            Console.WriteLine($"{0}\t {x:f1}\t {f(x):f5} \t{s:f5}");
            for (int i = 1; i < n - 1; i++)
            {
                x = x + h;
                s = s + f(x);
                Console.Write($"{i + 1}\t {x:f1}\t {f(x):f5} \t{s:f5}");
                Console.WriteLine($" \t {s * h:f5}");
            }
            I = s * h;
        }

        /// <summary>
        /// 变步长梯形求积公式(具有二阶精度)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="epsilon"></param>
        /// <param name="I"></param>
        public static void Integral_Trapezoidal_InvariableStep(F f, double a, double b
            , double epsilon //相对误差
            , out double I)
        {
            I = double.NaN;
            double h = b - a;
            double s1 = 0;
            double s2 = (f(a) + f(b)) * h / 2;
            double x = a;
            int n = 1;
            Console.WriteLine($"{n}\t {s2:f5}\t {(s2 - s1) / s2:f5}");
            while (Math.Abs((s2 - s1) / s2) > epsilon)
            {
                h = h / 2;
                s1 = s2;
                double T = 0;
                x = a + h;
                n = 2 * n;
                while (x < b)
                {
                    T = T + f(x);
                    x = x + 2 * h;
                }
                s2 = s1 / 2 + h * T;
                Console.WriteLine($"{n}\t {s2:f5}\t {(s2 - s1) / s2:f5}");
            }
            I = s2;
        }

        /// <summary>
        /// 辛普森求积公式 抛物线法(具有三阶精度)
        /// 1.将积分区间[a,b]分为2n个等间距的小区间
        /// 2.在每两个相邻小区间内，过曲线上的三个点做抛物线，抛物线的对称轴与Y轴平行
        /// 3.n个抛物线梯形的面积之和为函数的定积分
        /// I = h/3 * [ f(x0) - f(x2n) + SUM(4*f(x{2i-1}) + 2*f(x{2i})) ] i=1,...n
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="h"></param>
        /// <param name="I"></param>
        public static void Integral_Simpson(F f, double a, double b
            , double h //步长
            , out double I)
        {
            I = double.NaN;
            double n = Math.Floor((b - a) / (2 * h));
            double x0 = a;
            double x2n = b;
            double s0 = f(x0);
            double s2n = f(x2n);
            double sum1 = 0, sum2 = 0;
            for (int i = 0; i < n; i++)
            {
                double x1 = x0 + (2 * (i + 1) - 1) * h;
                double x2 = x1 + h;
                sum1 = sum1 + f(x1);
                sum2 = sum2 + f(x2);
            }
            double si = (s0 - s2n + 4 * sum1 + 2 * sum2) * h / 3;
            I = si;
        }

        /// <summary>
        /// 辛普森求积公式 抛物线法(具有三阶精度)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="h"></param>
        /// <param name="I"></param>
        public static void Integral_Simpson_Optimize(F f, double a, double b
            , double h //步长
            , out double I)
        {
            I = double.NaN;
            double n = Math.Floor((b - a) / h);
            double x0 = a;
            double xn = b;
            double s0 = f(x0);
            double sn = f(xn);
            double sum1 = 0, sum2 = 0;
            for (int i = 0; i < n; i += 2)
            {
                double x1 = x0 + (i + 1) * h;
                double x2 = x1 + h;
                sum1 = sum1 + f(x1);
                sum2 = sum2 + f(x2);
            }
            double si = (s0 - sn + 4 * sum1 + 2 * sum2) * h / 3;
            I = si;
        }

        /// <summary>
        /// 变步长辛普森求积公式
        /// h0 = (b-a) RP0 =f(a)+f(b) RC0 = SUM(f((a-h0/2)+i*h0)) i=1
        /// I1=h0/6 * (RP0 + 4 * RC0)
        /// hk = h{k-1}/2 = (b-a)/(2^k)  RPk = f(a) + f(b) + 2*SUM(f(a+i*hk))=RP{k-1} + 2RC{k-1} i=1,...,2n-1
        /// RCk = SUM(f((a-hk/2)+i*hk)) i=1,...,2n
        /// I2n = hk/6*(RPk+4*RCk)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="epsilon"></param>
        /// <param name="I"></param>
        public static void Integral_Simpson_InvariableStep(F f, double a, double b
            , double epsilon //相对误差
            , out double I)
        {
            I = double.NaN;
            int n = 1;
            double s1 = 0;
            double h = b - a;
            double rp = f(a) + f(b);
            double x = a - h / 2;
            double rc = f(x + h);
            double s2 = (rp + 4 * rc) * h / 6;
            Console.WriteLine($"{s2:f5}\t {n}\t {(s2 - s1) / s2:f5}");
            while (Math.Abs((s2 - s1) / s2) > epsilon)
            {
                s1 = s2;
                h = h / 2;
                n = n + n;
                rp = rp + 2 * rc;
                x = a - h / 2;
                rc = 0;
                for (int i = 0; i < n; i++)
                {
                    x = x + h;
                    rc = rc + f(x);
                }
                s2 = (rp + 4 * rc) * h / 6;
                Console.WriteLine($"{s2:f5}\t {n}\t {(s2 - s1) / s2:f5}");
            }
            I = s2;
        }

        /// <summary>
        /// 牛顿-柯特斯求积公式（具有四阶精度）
        /// 过曲线上的四个点做三次多项式曲线
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="h"></param>
        /// <param name="I"></param>
        public static void Integral_NewtonCotes(F f, double a, double b
            , double h
            , out double I)
        {
            I = double.NaN;
            int n = (int)Math.Floor((b - a) / h);
            double x0 = a;
            double xn = b;
            double s0 = f(x0);
            double sn = f(xn);
            double sum1 = 0, sum2 = 0, sum3 = -sn;
            for (int i = 0; i < n; i += 3)
            {
                double x1 = x0 + (i + 1) * h;
                double x2 = x1 + h;
                double x3 = x2 + h;
                sum1 += f(x1);
                sum2 += f(x2);
                sum3 += f(x3);
            }
            I = (s0 + sn + 3 * sum1 + 3 * sum2 + 2 * sum3) * 3 * h / 8;
        }

        public static void Integral_NewtonCotes_N(F f, double a, double b
            , double h
            , int _N
            , out double I)
        {
            I = double.NaN;
            int n = (int)Math.Floor((b - a) / h);
            double x0 = a;
            double xn = b;
            double s0 = f(x0);
            double sn = f(xn);
            var coes = FindNC_COEs(_N);
            double[] sums = new double[coes.Length];
            for (int i = 0; i < n; i += coes.Length)
            {
                double[] xs = new double[coes.Length];
                xs[0] = x0 + (i + 1) * h;
                for (int j = 1; j < xs.Length; j++)
                {
                    xs[j] = xs[j - 1] + h;
                }
                for (int j = 0; j < sums.Length; j++)
                {
                    sums[j] += f(xs[j]);
                }
            }
            I = 0;
            for (int i = 0; i < sums.Length; i++)
            {
                I += sums[i] * coes[i];
            }
            I = I * (b - a);
        }

        /// <summary>
        /// C_i^n 柯特斯系数表
        /// </summary>
        /// <param name="n"></param>
        /// <returns></returns>
        private static double[] FindNC_COEs(int n)
        {
            switch (n)
            {
                case 1:
                    return new double[] { 1d / 2d, 1d / 2d };
                case 2:
                    return new double[] { 1d / 6d, 2d / 3d, 1d / 6d };
                case 3:
                    return new double[] { 1d / 8d, 3d / 8d, 3d / 8d, 1d / 8d };
                case 4:
                    return new double[] { 7d / 90d, 16d / 45d, 2d / 15d, 16d / 45d, 7d / 90d };
                default:
                    return null;
            }
        }

        /// <summary>
        /// 龙贝格求积公式
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="eps"></param>
        /// <param name="I"></param>
        public static void Integral_Romberg(F f, double a, double b
            , double eps
            , out double I
            )
        {
            int N = 10;
            double[,] t = new double[N, N];
            I = double.NaN;
            double h = b - a;
            int k = 0;
            int m = 0;
            t[k, m] = (f(a) + f(b)) * h / 2d;
            double s1 = 0;
            double s2 = t[k, m];
            eps = 1e-6;
            Console.WriteLine($"k\t s1\t s2\t (s2-s1)/s2\t ");
            Console.WriteLine($"{k}\t {s1:f6}\t {s2:f6}\t {(s2 - s1) / s2:f7}\t ");
            while (Math.Abs((s2 - s1) / s2) > eps)
            {
                k++;
                if (k >= N) { break; }
                h = h / 2;
                double s = 0;
                for (int i = 1; i <= Math.Pow(2, k) - 1; i += 2)
                {
                    s = s + f(a + i * h);
                }
                t[k, m] = t[k - 1, m] / 2 + h * s;
                //龙贝格公式求：t[k: -1:1, 2:k + 1],并重新确定s1,s2
                for (int m1 = m; m1 < k; m1++)
                {
                    int k1 = k - m1 - 1;
                    double xx = (Math.Pow(4, m1 + 1) * t[k1 + 1, m1] - t[k1, m1])
                        / (Math.Pow(4, m1 + 1) - 1);
                    t[k1, m1 + 1] = xx;
                }
                s1 = s2;
                s2 = t[0, k];
                Console.WriteLine($"{k}\t {s1:f6}\t {s2:f6}\t {(s2 - s1) / s2:f7}\t ");
            }
            I = s2;
            Console.WriteLine();
            t.Print(6);
            Console.WriteLine();
            Console.WriteLine($"{k}\t {s1:f5}\t {s2:f5}\t {(s2 - s1) / s2:f5}\t ");
        }

        /// <summary>
        /// 方法一
        /// 在x方向将被积区间[a,b]划分为n等份，在y方向将被积区间[c,d]划分为m等份（n,m都是偶数）
        /// 由此将积分区域划分为 n x m 个小矩形，边长：h=(b-a)/n,k=(d-c)/m
        /// 各等距节点上的函数值：f(xi,yi)=fij=f(a+ih,c+jk)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="c"></param>
        /// <param name="d"></param>
        /// <param name="n"></param>
        /// <param name="m"></param>
        /// <param name="I"></param>
        public static void Integral_Simpson_xy(FF f
            , double a, double b
            , double c, double d
            , int n, int m
            , out double I
            )
        {
            I = double.NaN;
            if (n % 1 != 0)
                n++;
            if (m % 1 != 0)
                m++;
            double h = (b - a) / n;
            double hk = (d - c) / m;
            double sum = 0;
            double y = c;
            double j = 0;
            double L = 1;
            double p = 1;
            while (j <= m)
            {
                double x = a;
                for (int i = 0; i < n - 1; i += 2)
                {
                    sum = sum + p * (f(x, y) + 4 * f(x + h, y) + f(x + 2 * h, y));
                    x += 2 * h;
                }
                j++;
                y += hk;
                if (j < m)
                {
                    if (L == 1)
                    {
                        p = 4;
                        L++;
                    }
                    else
                    {
                        p = 2;
                        L = 1;
                    }
                }
                else
                {
                    p = 1.0;
                }
            }
            sum = sum * h * hk / 9;
            I = sum;
        }

        /// <summary>
        /// 方法二
        /// 对二重积分的区间[a,b],[c,d]进行与方法一相同的离散划分
        /// 两次利用辛普森求积公式，得到函数f(x,y)的二重数值积分
        /// </summary>
        /// <param name="f"></param>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="c"></param>
        /// <param name="d"></param>
        /// <param name="n"></param>
        /// <param name="m"></param>
        /// <param name="I"></param>
        public static void Integral_Simpson_xy2(FF f
            , double a, double b
            , double c, double d
            , int n, int m
            , out double I
            )
        {
            I = double.NaN;
            if (n % 1 != 0)
                n++;
            if (m % 1 != 0)
                m++;
            double[,] fxy = new double[n + 1, m + 1];
            double hk = (d - c) / m;
            int j = 0;
            double sum = 0;
            double si = 0;
            Simpsons(f, a, b, c, d, n, m, fxy, j, c, out si);
            sum += si;
            double y = 0;
            for (j = 1; j < m; j += 2)
            {
                y = c + j * hk;
                Simpsons(f, a, b, c, d, n, m, fxy, j, y, out si);
                sum += 4 * si;
                y = y + hk;
                Simpsons(f, a, b, c, d, n, m, fxy, j + 1, y, out si);
                sum += 2 * si;
            }
            j = m;
            Simpsons(f, a, b, c, d, n, m, fxy, j, d, out si);
            sum = (sum - si) * hk / 3;
            fxy.Print();
            I = sum;
        }

        private static void Simpsons(FF f
            , double a, double b, double c, double d
            , int n, int m, double[,] fxy
            , int j
            , double y
            , out double si
            )
        {
            si = double.NaN;
            int i = 0;
            double x0 = a;
            double h = (b - a) / n;
            double s0 = f(x0, y);
            fxy[i, j] = s0;
            double xn = b;
            double sn = f(xn, y);
            fxy[n, j] = sn;
            double sum1 = 0;
            double sum2 = -sn;
            for (i = 1; i <= n; i += 2)
            {
                double x1 = x0 + i * h;
                double x2 = x1 + h;
                fxy[i, j] = f(x1, y);
                sum1 += fxy[i, j];
                fxy[i + 1, j] = f(x2, y);
                sum2 += fxy[i + 1, j];
            }
            si = (s0 + sn + 4 * sum1 + 2 * sum2) * h / 3;
        }

        /// <summary>
        /// 利用插值多项式求数值微分
        /// 1.两点公式
        /// f'(x)=[f(x1)-f(x2)]/(x1-x2)
        /// 2.三点公式
        /// f'(x)=[-3f(x1)+4f(x2)-f(x3)]/(2h)
        /// f''(x)=[f(x0)-3f(x1)+f(x2)]/(h^2)
        /// </summary>
        /// <param name="f"></param>
        /// <param name="x"></param>
        /// <param name="h"></param>
        /// <param name="eps"></param>
        /// <param name="dx"></param>
        /// <param name="dxx"></param>
        public static void Difference_Finite(F f
            , double x
            , double h
            , double eps
            , out double dx
            , out double dxx)
        {
            double f0 = f(x);
            double dx2 = double.NaN, dx3 = double.NaN, dxx3 = double.NaN;
            while (h >= eps 
                && h > 1e-7) //计算误差
            {
                double f1 = f(x + h);
                double f2 = f(x + 2.0 * h);
                dx2 = (f1 - f0) / h;
                dx3 = (-3.0 * f0 + 4.0 * f1 - f2) / 2.0 / h;
                var da = (f2 + f0 - 2.0 * f1);
                var db = h * h;
                dxx3 = (f2 - 2.0 * f1 + f0) / (h * h);
                var dxx4 = (f2 / db - 2.0 * f1 / db + f0 / db);
                h = h / 10.0;
                //Console.WriteLine($"{f0} \t {f1} \t {f2} \t {db} \t {dxx3}");
                Console.WriteLine($"{h:f6}\t {dx2:f6}\t {dx3:f6}\t {dxx3:f6}");
            }
            dx = dx3;
            dxx = dxx3;
        }
    }
}
